/**
 * 
 */
package ua.lviv.franko.solvers.stressanalysis.oneDimention;

import java.util.ArrayList;
import java.util.Iterator;
import java.util.List;

import ua.lviv.franko.integration.ConstantFunction1D;
import ua.lviv.franko.integration.ConstantFunction2D;
import ua.lviv.franko.integration.GaussQuadrature;
import ua.lviv.franko.integration.IFunction1D;
import ua.lviv.franko.integration.IFunction2D;
import ua.lviv.franko.integration.Quadrature;
import ua.lviv.franko.matrix.RowType;
import ua.lviv.franko.matrix.hash.gauss.Gauss;
import ua.lviv.franko.matrix.hash.gauss.HashMatrixGaussStyle;
import ua.lviv.franko.matrix.hash.lu.HashMatrixLUStyle;
import ua.lviv.franko.matrix.hash.lu.LUFactorization;
import ua.lviv.franko.mlsFunctions.oneDimension.BaseFunction;
import ua.lviv.franko.mlsFunctions.oneDimension.Particle;
import ua.lviv.franko.mlsFunctions.oneDimension.Values;
import ua.lviv.franko.mlsFunctions.oneDimension.weightFunctions.CubicSpline;
import ua.lviv.franko.mlsFunctions.oneDimension.weightFunctions.IWeightFunction;
import ua.lviv.franko.visualisation.EFGResult1D;

/**
 * @author marko.bekhta
 * 
 */
public class TimoshenkoMembraneMFSolver {

	// general data
	protected static final int		QUADRATURE_POINTS_NUMBER		= 7;
	protected static final int		NUMBER_OF_INTEGRATION_INTERVALS	= 4;
	protected static final int		NUMBER_OF_UNKNOWNS				= 3;
	protected Quadrature			gau;
	protected EquationParameters	parameters;
	protected HashMatrixLUStyle		mainMatrix;
	protected double[]				b;
	protected double				leftBoundery;
	protected double				rightBoundery;
	protected BoundaryCondition		leftBC;
	protected BoundaryCondition		rightBC;

	// problem data
	protected double				kKk1;
	protected double				kKk2;
	protected double				kKk3;
	protected IFunction1D			rho1;
	protected IFunction1D			rho3;
	protected IFunction1D			mu1;

	// data for base functions
	protected List<Particle>		particles;
	protected List<BaseFunction>	baseFunctions;

	/**
	 * 
	 * @param parameters
	 *            - equation parameters
	 * @param particles
	 *            - list of particles
	 * @param leftBoundery
	 *            - left side of domain
	 * @param rightBoundery
	 *            - right side of domain
	 * @param PolinomialDegree
	 *            - degree of polynomials used for basis
	 * @param weight
	 *            - weight function
	 */
	public TimoshenkoMembraneMFSolver(final EquationParameters parameters, List<Particle> particles, double leftBoundery, double rightBoundery, int PolinomialDegree,
			IWeightFunction weight) {

		super();
		this.parameters = parameters;
		this.particles = particles;
		this.rightBoundery = rightBoundery;
		this.leftBoundery = leftBoundery;

		gau = new GaussQuadrature(QUADRATURE_POINTS_NUMBER);

		kKk1 = parameters.getE() * parameters.getH() / (1 - parameters.getV() * parameters.getV());
		kKk2 = parameters.getkPrime() * parameters.getG() * parameters.getH();
		kKk3 = parameters.getE() * Math.pow(parameters.getH(), 3) / (12 * (1 - parameters.getV() * parameters.getV()));
		rho1 = new IFunction1D() {

			@Override
			public double calculate(final double xxx) {
				double h = parameters.getH();
				final double k1 = parameters.getK1().calculate(xxx);

				return gau.integrate(-h / 2, h / 2, new IFunction1D() {

					@Override
					public double calculate(double x) {
						double f1 = parameters.getF1().calculate(xxx, x);
						return (1 + k1 * x) * f1;
					}

					@Override
					public double calculateDerivaty(double x) {

						return 0;
					}
				});
			}

			@Override
			public double calculateDerivaty(double xxx) {

				return 0;
			}
		};

		rho3 = new IFunction1D() {

			@Override
			public double calculate(final double xxx) {
				double h = parameters.getH();
				final double k1 = parameters.getK1().calculate(xxx);

				return gau.integrate(-h / 2, h / 2, new IFunction1D() {

					@Override
					public double calculate(double x) {
						double f1 = parameters.getF3().calculate(xxx, x);
						return (1 + k1 * x) * f1;
					}

					@Override
					public double calculateDerivaty(double x) {

						return 0;
					}
				});
			}

			@Override
			public double calculateDerivaty(double xxx) {

				return 0;
			}
		};

		mu1 = new IFunction1D() {

			@Override
			public double calculate(final double xxx) {
				double h = parameters.getH();
				final double k1 = parameters.getK1().calculate(xxx);

				return gau.integrate(-h / 2, h / 2, new IFunction1D() {

					@Override
					public double calculate(double x) {
						double f1 = parameters.getF1().calculate(xxx, x);
						return (1 + k1 * x) * x * f1;
					}

					@Override
					public double calculateDerivaty(double x) {

						return 0;
					}
				});
			}

			@Override
			public double calculateDerivaty(double xxx) {

				return 0;
			}
		};
		// creation of base functions;
		this.createBaseFunctions(PolinomialDegree, weight);
	}

	public void generateMatrix() {
		int tmpN = 0;
		if (leftBC.isDirihlet()) {
			tmpN++;
		}
		if (rightBC.isDirihlet()) {
			tmpN++;
		}
		mainMatrix = new HashMatrixLUStyle(RowType.ARRAY_TYPE, (particles.size() + tmpN) * NUMBER_OF_UNKNOWNS);
		b = new double[(particles.size() + tmpN) * NUMBER_OF_UNKNOWNS];

		for (int i = 0; i < particles.size(); i++) {
			// System.out.println(i);
			BaseFunction baseFunction = baseFunctions.get(i);
			double[][] vals = this.gaussQuadratureIntegration(NUMBER_OF_INTEGRATION_INTERVALS,
					Math.max(baseFunction.getXi().getXi() - baseFunction.getXi().getRho(), leftBoundery),
					Math.min(baseFunction.getXi().getXi() + baseFunction.getXi().getRho(), rightBoundery), baseFunction);
			int pos1 = particles.indexOf(baseFunction.getXi());
			for (int j = 0; j < baseFunction.getNodeList().size(); j++) {
				int pos2 = particles.indexOf(baseFunction.getNodeList().get(j));
				mainMatrix.setValue(pos1 * 3, pos2 * 3, vals[j * NUMBER_OF_UNKNOWNS][0]);
				mainMatrix.setValue(pos1 * 3, pos2 * 3 + 1, vals[j * NUMBER_OF_UNKNOWNS][1]);
				mainMatrix.setValue(pos1 * 3, pos2 * 3 + 2, vals[j * NUMBER_OF_UNKNOWNS][2]);

				mainMatrix.setValue(pos1 * 3 + 1, pos2 * 3, vals[j * NUMBER_OF_UNKNOWNS + 1][0]);
				mainMatrix.setValue(pos1 * 3 + 1, pos2 * 3 + 1, vals[j * NUMBER_OF_UNKNOWNS + 1][1]);
				mainMatrix.setValue(pos1 * 3 + 1, pos2 * 3 + 2, vals[j * NUMBER_OF_UNKNOWNS + 1][2]);

				mainMatrix.setValue(pos1 * 3 + 2, pos2 * 3, vals[j * NUMBER_OF_UNKNOWNS + 2][0]);
				mainMatrix.setValue(pos1 * 3 + 2, pos2 * 3 + 1, vals[j * NUMBER_OF_UNKNOWNS + 2][1]);
				mainMatrix.setValue(pos1 * 3 + 2, pos2 * 3 + 2, vals[j * NUMBER_OF_UNKNOWNS + 2][2]);

				mainMatrix.setValue(pos2 * 3, pos1 * 3, vals[j * NUMBER_OF_UNKNOWNS][3]);
				mainMatrix.setValue(pos2 * 3, pos1 * 3 + 1, vals[j * NUMBER_OF_UNKNOWNS][4]);
				mainMatrix.setValue(pos2 * 3, pos1 * 3 + 2, vals[j * NUMBER_OF_UNKNOWNS][5]);

				mainMatrix.setValue(pos2 * 3 + 1, pos1 * 3, vals[j * NUMBER_OF_UNKNOWNS + 1][3]);
				mainMatrix.setValue(pos2 * 3 + 1, pos1 * 3 + 1, vals[j * NUMBER_OF_UNKNOWNS + 1][4]);
				mainMatrix.setValue(pos2 * 3 + 1, pos1 * 3 + 2, vals[j * NUMBER_OF_UNKNOWNS + 1][5]);

				mainMatrix.setValue(pos2 * 3 + 2, pos1 * 3, vals[j * NUMBER_OF_UNKNOWNS + 2][3]);
				mainMatrix.setValue(pos2 * 3 + 2, pos1 * 3 + 1, vals[j * NUMBER_OF_UNKNOWNS + 2][4]);
				mainMatrix.setValue(pos2 * 3 + 2, pos1 * 3 + 2, vals[j * NUMBER_OF_UNKNOWNS + 2][5]);

			}
		}
	}

	public BoundaryCondition getLeftBC() {
		return leftBC;
	}

	public BoundaryCondition getRightBC() {
		return rightBC;
	}

	public void setLeftBC(BoundaryCondition leftBC) {
		this.leftBC = leftBC;
	}

	public void setRightBC(BoundaryCondition rightBC) {
		this.rightBC = rightBC;
	}

	public EFGResult1D solve() {
		if (leftBC.isDirihlet()) {
			BaseFunction bf = baseFunctions.get(0);
			// TODO
			double[] vals = bf.calculateAllFunctionsValuesOnly(leftBoundery);
			// double[] vals = bf.calculateAllFunctionsValuesOnlyGauss(leftBoundery);
			int size = mainMatrix.getN();
			for (int i = 0; i < vals.length; i++) {
				mainMatrix.setValue(size - 3, i * NUMBER_OF_UNKNOWNS, vals[i]);
				mainMatrix.setValue(size - 2, i * NUMBER_OF_UNKNOWNS + 1, vals[i]);
				mainMatrix.setValue(size - 1, i * NUMBER_OF_UNKNOWNS + 2, vals[i]);

				mainMatrix.setValue(i * NUMBER_OF_UNKNOWNS, size - 3, vals[i]);
				mainMatrix.setValue(i * NUMBER_OF_UNKNOWNS + 1, size - 2, vals[i]);
				mainMatrix.setValue(i * NUMBER_OF_UNKNOWNS + 2, size - 1, vals[i]);
			}

			b[size - 3] = leftBC.uVal;
			b[size - 2] = leftBC.vVal;
			b[size - 1] = leftBC.mVAl;

			if (rightBC.isDirihlet()) {
				bf = baseFunctions.get(baseFunctions.size() - 1);
				vals = bf.calculateAllFunctionsValuesOnly(rightBoundery);
				for (int i = 0; i < vals.length; i++) {
					mainMatrix.setValue(size - 6, (i + baseFunctions.size() - vals.length) * NUMBER_OF_UNKNOWNS, vals[i]);
					mainMatrix.setValue(size - 5, (i + baseFunctions.size() - vals.length) * NUMBER_OF_UNKNOWNS + 1, vals[i]);
					mainMatrix.setValue(size - 4, (i + baseFunctions.size() - vals.length) * NUMBER_OF_UNKNOWNS + 2, vals[i]);

					mainMatrix.setValue((i + baseFunctions.size() - vals.length) * NUMBER_OF_UNKNOWNS, size - 6, vals[i]);
					mainMatrix.setValue((i + baseFunctions.size() - vals.length) * NUMBER_OF_UNKNOWNS + 1, size - 5, vals[i]);
					mainMatrix.setValue((i + baseFunctions.size() - vals.length) * NUMBER_OF_UNKNOWNS + 2, size - 4, vals[i]);
				}
				b[size - 6] = rightBC.uVal;
				b[size - 5] = rightBC.vVal;
				b[size - 4] = rightBC.mVAl;
			}
		} else {
			if (rightBC.isDirihlet()) {
				BaseFunction bf = baseFunctions.get(baseFunctions.size() - 1);
				double[] vals = bf.calculateAllFunctionsValuesOnly(rightBoundery);
				int size = mainMatrix.getN();
				for (int i = 0; i < vals.length; i++) {
					mainMatrix.setValue(size - 6, (i + baseFunctions.size() - vals.length) * NUMBER_OF_UNKNOWNS, vals[i]);
					mainMatrix.setValue(size - 5, (i + baseFunctions.size() - vals.length) * NUMBER_OF_UNKNOWNS + 1, vals[i]);
					mainMatrix.setValue(size - 4, (i + baseFunctions.size() - vals.length) * NUMBER_OF_UNKNOWNS + 2, vals[i]);
				}
				b[size - 6] = rightBC.uVal;
				b[size - 5] = rightBC.vVal;
				b[size - 4] = rightBC.mVAl;

			}
		}

		// double val = 1.0;
		// b[0] = val;
		// b[1] = val;
		// b[2] = val;
		// val = 0.0;
		// b[size - 3] = val;
		// b[size - 2] = val;
		// b[size - 1] = val;

		// for (int i = 0; i < mainMatrix.getN(); i++) {
		// for (int j = 0; j < mainMatrix.getN(); j++) {
		// System.out.print(mainMatrix.getValue(i, j) + "\t");
		// }
		// System.out.println();// System.out.println(b[i]);
		// }
		// System.out.println("-------------\n-------------");
		HashMatrixLUStyle.printMatrix(mainMatrix);
		LUFactorization lu = new LUFactorization(this.mainMatrix);
		double[] results = lu.calculate(b);
		// for (int i = 0; i < particles.size(); i++) {
		// System.out.println(results[i * NUMBER_OF_UNKNOWNS]);
		// }
		// System.out.println("-------------\n-------------");
		return new EFGResult1D(results, baseFunctions, particles);
	}

	private void createBaseFunctions(int polinomialDegree, IWeightFunction weight) {
		baseFunctions = new ArrayList<BaseFunction>();
		for (Iterator<Particle> iterator = particles.iterator(); iterator.hasNext();) {
			Particle particle = iterator.next();
			BaseFunction bf = new BaseFunction(polinomialDegree, particle, weight);
			bf.setNodeList(this.getParticleSupport(particle, weight));
			this.baseFunctions.add(bf);
		}

	}

	private double[][] gaussQuadratureIntegration(int n, double ax, double bx, BaseFunction baseFunction) {
		double hx = (bx - ax) / n;

		double[][] result = new double[baseFunction.getNodeList().size() * NUMBER_OF_UNKNOWNS][2 * NUMBER_OF_UNKNOWNS];

		for (int i = 0; i < n; i++) {
			double[][] tmp = gaussQuadratureIntegration(ax + i * hx, ax + (i + 1) * hx, baseFunction);
			for (int k = 0; k < tmp.length; k++) {
				for (int q = 0; q < tmp[0].length; q++) {
					result[k][q] += tmp[k][q];
				}
			}
		}

		return result;
	}

	private double[][] gaussQuadratureIntegration(double ax, double bx, BaseFunction baseFunction) {
		// first column - baseFunction*v; second column - v*baseFunction
		double[][] result = new double[baseFunction.getNodeList().size() * NUMBER_OF_UNKNOWNS][2 * NUMBER_OF_UNKNOWNS];
		// its a position of main node - node of base function
		int pos = baseFunction.getNodeList().indexOf(baseFunction.getXi());

		double[] f = new double[NUMBER_OF_UNKNOWNS];

		for (int i = 0; i < QUADRATURE_POINTS_NUMBER; i++) {
			double xi = (ax + bx) / 2 + (bx - ax) / 2 * gau.getX(i);
			// TODO : change this ....
			Values[] values = baseFunction.calculateAllFunctions(xi);
			// Values[] values = baseFunction.calculateAllFunctionsGauss(xi);
			double fiDer = values[pos].valDerivaty;
			double fiVal = values[pos].val;
			double a1 = parameters.getA1().calculate(xi);
			double k1 = parameters.getK1().calculate(xi);

			double h = parameters.getH() / 2.0;

			double rho1Val = rho1.calculate(xi);
			double rho3Val = rho3.calculate(xi);
			double mu1Val = mu1.calculate(xi);

			double sigma11Minus = parameters.getSigma11Minus().calculate(xi);
			double sigma11Plus = parameters.getSigma11Plus().calculate(xi);

			double sigma13Minus = parameters.getSigma13Minus().calculate(xi);
			double sigma13Plus = parameters.getSigma13Plus().calculate(xi);

			double sigma33Minus = parameters.getSigma33Minus().calculate(xi);
			double sigma33Plus = parameters.getSigma33Plus().calculate(xi);

			for (int j = 0; j < values.length; j++) {
				double fjDer = values[j].valDerivaty;
				double fjVal = values[j].val;

				double[] res1 = new double[] {
						// first row
						kKk1 / a1 * fiDer * fjDer + a1 * k1 * k1 * kKk2 * fiVal * fjVal, kKk1 / a1 * k1 * fjVal * fiDer - k1 * kKk2 * fjDer * fiVal,
						-a1 * k1 * kKk2 * fiVal * fjVal,
						// second row
						-kKk2 * k1 * fjVal * fiDer + k1 * kKk1 * fjDer * fiVal, kKk2 / a1 * fjDer * fiDer + a1 * k1 * k1 * kKk1 * fiVal * fjVal, kKk2 * fjVal * fiDer,
						// third row
						-a1 * kKk2 * k1 * fiVal * fjVal, kKk2 * fjDer * fiVal, kKk3 / a1 * fjDer * fiDer + a1 * kKk2 * fiVal * fjVal };

				double[] res2 = new double[] {
						// first row
						kKk1 / a1 * fjDer * fiDer + a1 * k1 * k1 * kKk2 * fjVal * fiVal, kKk1 / a1 * k1 * fiVal * fjDer - k1 * kKk2 * fiDer * fjVal,
						-a1 * k1 * kKk2 * fjVal * fiVal,
						// second row
						-kKk2 * k1 * fiVal * fjDer + k1 * kKk1 * fiDer * fjVal, kKk2 / a1 * fiDer * fjDer + a1 * k1 * k1 * kKk1 * fjVal * fiVal, kKk2 * fiVal * fjDer,
						// thjrd row
						-a1 * kKk2 * k1 * fjVal * fiVal, kKk2 * fiDer * fjVal, kKk3 / a1 * fiDer * fjDer + a1 * kKk2 * fjVal * fiVal };
				for (int k = 0; k < NUMBER_OF_UNKNOWNS; k++) {
					result[j * 3 + k][0] += gau.getW(i) * res1[NUMBER_OF_UNKNOWNS * k];
					result[j * 3 + k][1] += gau.getW(i) * res1[NUMBER_OF_UNKNOWNS * k + 1];
					result[j * 3 + k][2] += gau.getW(i) * res1[NUMBER_OF_UNKNOWNS * k + 2];

					result[j * 3 + k][3] += gau.getW(i) * res2[NUMBER_OF_UNKNOWNS * k];
					result[j * 3 + k][4] += gau.getW(i) * res2[NUMBER_OF_UNKNOWNS * k + 1];
					result[j * 3 + k][5] += gau.getW(i) * res2[NUMBER_OF_UNKNOWNS * k + 2];
				}
			}
			f[0] += gau.getW(i) * fiVal * ((1 + k1 * h) * sigma13Plus + (1 - k1 * h) * sigma13Minus + rho1Val);
			f[1] += gau.getW(i) * fiVal * ((1 + k1 * h) * sigma33Plus + (1 - k1 * h) * sigma33Minus + rho3Val);
			f[2] += gau.getW(i) * fiVal * (((1 + k1 * h) * sigma11Plus - (1 - k1 * h) * sigma11Minus) * h + mu1Val);

			// f += gau.getW(i) * values[pos].val;
		}
		double bMinusA = (bx - ax) / 2.0;
		// f *= (bx - ax) / 2.0;
		f[0] *= bMinusA;
		f[1] *= bMinusA;
		f[2] *= bMinusA;

		for (int i = 0; i < result.length; i++) {
			result[i][0] *= bMinusA;
			result[i][1] *= bMinusA;
			result[i][2] *= bMinusA;
			result[i][3] *= bMinusA;
			result[i][4] *= bMinusA;
			result[i][5] *= bMinusA;
		}

		b[particles.indexOf(baseFunction.getXi()) * NUMBER_OF_UNKNOWNS] = f[0];
		b[particles.indexOf(baseFunction.getXi()) * NUMBER_OF_UNKNOWNS + 1] = f[1];
		b[particles.indexOf(baseFunction.getXi()) * NUMBER_OF_UNKNOWNS + 2] = f[2];
		return result;
	}

	private List<Particle> getParticleSupport(Particle particle, IWeightFunction weight) {
		List<Particle> list = new ArrayList<Particle>();
		for (Iterator<Particle> iterator = particles.iterator(); iterator.hasNext();) {
			Particle p = (Particle) iterator.next();
			if (weight.intersectedSupport(particle, p)) {
				list.add(p);
			}
		}
		return list;
		// return particles;
	}

	/**
	 * @param args
	 */
	public static void main(String[] args) {

		IFunction2D f1 = new ConstantFunction2D(1.0);
		IFunction2D f3 = new ConstantFunction2D(1.0);
		IFunction1D k1 = new ConstantFunction1D(0.0);
		IFunction1D a1 = new ConstantFunction1D(1.0);
		IFunction1D sigma11Plus = new ConstantFunction1D(0.0);
		IFunction1D sigma13Minus = new ConstantFunction1D(0.0);
		IFunction1D sigma11Minus = new ConstantFunction1D(0.0);
		IFunction1D sigma33Minus = new ConstantFunction1D(0.0);
		IFunction1D sigma33Plus = new ConstantFunction1D(0.0);
		IFunction1D sigma13Plus = new ConstantFunction1D(0.0);
		double v = 0.3;
		double e = 10000;
		double h = 0.5;

		EquationParameters eqP = new EquationParameters(a1, k1, e, v, h, f1, f3, sigma13Plus, sigma13Minus, sigma11Plus, sigma11Minus, sigma33Plus, sigma33Minus);
		double a = 0.0;
		double b = 1.0;
		int n = 10;
		double hStep = (b - a) / (double) n;
		double rho = 1.50000001 * hStep;
		List<Particle> particles = new ArrayList<Particle>();

		for (int i = 0; i <= n; i++) {
			particles.add(new Particle(a + i * hStep, rho));
		}

		TimoshenkoMembraneMFSolver s = new TimoshenkoMembraneMFSolver(eqP, particles, 0.0, 1.0, 1, new CubicSpline());
		s.setLeftBC(new BoundaryCondition(0.0, 0.0, 0.0, true));
		s.setRightBC(new BoundaryCondition(0.0, 0.0, 0.0, true));
		long start = System.nanoTime();
		s.generateMatrix();
		// System.err.println("generateMatrix() time is = " + (System.nanoTime() - start));
		// start = System.nanoTime();
		EFGResult1D res = s.solve();
		System.err.println("solve() time is = " + (System.nanoTime() - start));
		res.display(a, b);
	}

}
